# PoP - Policing Socio-Geographic Boundaries and Inequality
# script for creating figures and tables in Appendix E
# Appendix E Table 11: Influence of Percent White on Logged Misdemeanor Arrests (Standardized), Conditional on Racial Boundary Status
# Appendix E Table 12: Influence of Percent White on Logged Felony Arrests (Standardized), Conditional on Racial Boundary Status
# Appendix E Figure 7: Influence of Percent Poverty on Logged Arrests by Type (Standardized), Condi- tional on Racial Boundary Status.
# Appendix E Table 13: Influence of Percent Poverty on Logged Misdemeanor Arrests (Standardized), Condi- tional on Racial Boundary Status
# Appendix E Table 14: Influence of Percent Poverty on Logged Felony Arrests (Standardized), Conditional on Racial Boundary Status
# Appendix E Figure 8:  Louisville: Influence of Percent White on Logged Arrests by Type (Standardized), Conditional on Racial Boundary Status.
# Appendix E Table 15:  Louisville: Influence of Percent White on Logged Arrests by Type (Standardized), Conditional on Racial Boundary Status.



suppressPackageStartupMessages(
  
  {
    library(AER)
    library(dplyr)
    library(fixest)
    library(lfe)
    library(tidyverse)
    library(ggplot2)
    library(MASS)
    library(sensemakr)
    library(haven)
    library(readstata13)
    library(readxl)
    library(readr)
    library(gridExtra)
    library(areal)
    library(car)
    library(estimatr)
    library(magrittr)
    library(texreg)
    library(sandwich)
    library(jtools)
    library(ggthemes)
    library(meta)
  }
)

# first prepare variables by city
# read in data by city

# Load Atlanta final dataset
load("atl_final.RData")

# log and +1 to variables (pop already logged)
atl_fin = atl_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
atl_fin$larrest_sd <- atl_fin$larrests - (mean(atl_fin$larrests)/sd(atl_fin$larrests))
atl_fin$lmisdemeanors_sd <- atl_fin$lmisdemeanors - (mean(atl_fin$lmisdemeanors)/sd(atl_fin$lmisdemeanors))
atl_fin$lfelonies_sd <- atl_fin$lfelonies - (mean(atl_fin$lfelonies)/sd(atl_fin$lfelonies))
atl_fin$lnonviolent_sd <- atl_fin$lnonviolent - (mean(atl_fin$lnonviolent)/sd(atl_fin$lnonviolent))
atl_fin$lviolent_sd <- atl_fin$lviolent - (mean(atl_fin$lviolent)/sd(atl_fin$lviolent))
atl_fin$lsociety_sd <- atl_fin$lsociety - (mean(atl_fin$lsociety)/sd(atl_fin$lsociety))
atl_fin$lperson_sd <- atl_fin$lperson - (mean(atl_fin$lperson)/sd(atl_fin$lperson))
atl_fin$lproperty_sd <- atl_fin$lproperty - (mean(atl_fin$lproperty)/sd(atl_fin$lproperty))

# rename racial and economic boundary variable
atl_fin$atl_white_blv <- atl_fin$p_race_white_blv
atl_fin$atl_ses_blv <- atl_fin$ses_blv

# Load Austin final dataset
load("aus_final.RData")



# log and +1 to variables (pop already logged)
aus_fin = aus_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
aus_fin$larrest_sd <- aus_fin$larrests - (mean(aus_fin$larrests)/sd(aus_fin$larrests))
aus_fin$lmisdemeanors_sd <- aus_fin$lmisdemeanors - (mean(aus_fin$lmisdemeanors)/sd(aus_fin$lmisdemeanors))
aus_fin$lfelonies_sd <- aus_fin$lfelonies - (mean(aus_fin$lfelonies)/sd(aus_fin$lfelonies))
aus_fin$lnonviolent_sd <- aus_fin$lnonviolent - (mean(aus_fin$lnonviolent)/sd(aus_fin$lnonviolent))
aus_fin$lviolent_sd <- aus_fin$lviolent - (mean(aus_fin$lviolent)/sd(aus_fin$lviolent))
aus_fin$lsociety_sd <- aus_fin$lsociety - (mean(aus_fin$lsociety)/sd(aus_fin$lsociety))
aus_fin$lperson_sd <- aus_fin$lperson - (mean(aus_fin$lperson)/sd(aus_fin$lperson))
aus_fin$lproperty_sd <- aus_fin$lproperty - (mean(aus_fin$lproperty)/sd(aus_fin$lproperty))

# rename racial and economic boundary variable
aus_fin$aus_white_blv <- aus_fin$p_race_white_blv
aus_fin$aus_ses_blv <- aus_fin$ses_blv


## Load Boston data
load("bos_final.RData")


# log and +1 to variables (pop already logged)
bos_fin = bos_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime + 1))



# create scaled DVs (to account for different pop. sizes and density in blocks)
bos_fin$larrests_sd <- bos_fin$larrests - (mean(bos_fin$larrests)/sd(bos_fin$larrests))
bos_fin$lmisdemeanors_sd <- bos_fin$lmisdemeanors - (mean(bos_fin$lmisdemeanors)/sd(bos_fin$lmisdemeanors))
bos_fin$lfelonies_sd <- bos_fin$lfelonies - (mean(bos_fin$lfelonies)/sd(bos_fin$lfelonies))
bos_fin$lnonviolent_sd <- bos_fin$lnonviolent - (mean(bos_fin$lnonviolent)/sd(bos_fin$lnonviolent))
bos_fin$lviolent_sd <- bos_fin$lviolent - (mean(bos_fin$lviolent)/sd(bos_fin$lviolent))
bos_fin$lsociety_sd <- bos_fin$lsociety - (mean(bos_fin$lsociety)/sd(bos_fin$lsociety))
bos_fin$lperson_sd <- bos_fin$lperson - (mean(bos_fin$lperson)/sd(bos_fin$lperson))
bos_fin$lproperty_sd <- bos_fin$lproperty - (mean(bos_fin$lproperty)/sd(bos_fin$lproperty))

# rename racial and economic boundary variable
bos_fin$bos_white_blv <- bos_fin$p_race_white_blv
bos_fin$bos_ses_blv <- bos_fin$ses_blv

## Load Chicago data
load("chi_final.RData")

# log and +1 to variables (pop already logged)
chi_fin = chi_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(property_crime+1),
         lviolentcrime = log(violent_crime+1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
chi_fin$larrest_sd <- chi_fin$larrests - (mean(chi_fin$larrests)/sd(chi_fin$larrests))
chi_fin$lmisdemeanors_sd <- chi_fin$lmisdemeanors - (mean(chi_fin$lmisdemeanors)/sd(chi_fin$lmisdemeanors))
chi_fin$lnonviolent_sd <- chi_fin$lnonviolent - (mean(chi_fin$lnonviolent)/sd(chi_fin$lnonviolent))
chi_fin$lsociety_sd <- chi_fin$lsociety - (mean(chi_fin$lsociety)/sd(chi_fin$lsociety))
chi_fin$lfelonies_sd <- chi_fin$lfelonies - (mean(chi_fin$lfelonies)/sd(chi_fin$lfelonies))
chi_fin$lperson_sd <- chi_fin$lperson - (mean(chi_fin$lperson)/sd(chi_fin$lperson))
chi_fin$lproperty_sd <- chi_fin$lproperty - (mean(chi_fin$lproperty)/sd(chi_fin$lproperty))
chi_fin$lviolent_sd <- chi_fin$lviolent - (mean(chi_fin$lviolent)/sd(chi_fin$lviolent))


# rename racial and economic boundary variable
chi_fin$chi_white_blv <- chi_fin$p_race_white_blv
chi_fin$chi_ses_blv <- chi_fin$ses_blv

## Load Louisville data
load("lou_final.RData")


# log and +1 to variables (pop already logged) (no misdemeanors vs felonies for louisville)
lou_fin = lou_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property+1),
         lviolentcrime = log(crime_violent+1),
         lviolent = log(violent_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
lou_fin$larrest_sd <- lou_fin$larrests - (mean(lou_fin$larrests)/sd(lou_fin$larrests))
lou_fin$lnonviolent_sd <- lou_fin$lnonviolent - (mean(lou_fin$lnonviolent)/sd(lou_fin$lnonviolent))
lou_fin$lsociety_sd <- lou_fin$lsociety - (mean(lou_fin$lsociety)/sd(lou_fin$lsociety))
lou_fin$lviolent_sd <- lou_fin$lviolent - (mean(lou_fin$lviolent)/sd(lou_fin$lviolent))
lou_fin$lperson_sd <- lou_fin$lperson - (mean(lou_fin$lperson)/sd(lou_fin$lperson))
lou_fin$lproperty_sd <- lou_fin$lproperty - (mean(lou_fin$lproperty)/sd(lou_fin$lproperty))

# rename racial and economic boundary variable
lou_fin$lou_white_blv <- lou_fin$p_race_white_blv
lou_fin$lou_ses_blv <- lou_fin$ses_blv

## Load Milwaukee data
load("mil_final.RData")


# log and +1 to variables (pop already logged)
mil_fin = mil_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
mil_fin$larrest_sd <- mil_fin$larrests - (mean(mil_fin$larrests)/sd(mil_fin$larrests))
mil_fin$lmisdemeanors_sd <- mil_fin$lmisdemeanors - (mean(mil_fin$lmisdemeanors)/sd(mil_fin$lmisdemeanors))
mil_fin$lfelonies_sd <- mil_fin$lfelonies - (mean(mil_fin$lfelonies)/sd(mil_fin$lfelonies))
mil_fin$lnonviolent_sd <- mil_fin$lnonviolent - (mean(mil_fin$lnonviolent)/sd(mil_fin$lnonviolent))
mil_fin$lviolent_sd <- mil_fin$lviolent - (mean(mil_fin$lviolent)/sd(mil_fin$lviolent))
mil_fin$lsociety_sd <- mil_fin$lsociety - (mean(mil_fin$lsociety)/sd(mil_fin$lsociety))
mil_fin$lperson_sd <- mil_fin$lperson - (mean(mil_fin$lperson)/sd(mil_fin$lperson))
mil_fin$lproperty_sd <- mil_fin$lproperty - (mean(mil_fin$lproperty)/sd(mil_fin$lproperty))

# rename racial and economic boundary variable
mil_fin$mil_white_blv <- mil_fin$p_race_white_blv
mil_fin$mil_ses_blv <- mil_fin$ses_blv


## Load Seattle data
load("sea_final.RData")

# log and +1 to variables (pop already logged)
sea_fin = sea_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
sea_fin$larrest_sd <- sea_fin$larrests - (mean(sea_fin$larrests)/sd(sea_fin$larrests))
sea_fin$lmisdemeanors_sd <- sea_fin$lmisdemeanors - (mean(sea_fin$lmisdemeanors)/sd(sea_fin$lmisdemeanors))
sea_fin$lfelonies_sd <- sea_fin$lfelonies - (mean(sea_fin$lfelonies)/sd(sea_fin$lfelonies))
sea_fin$lnonviolent_sd <- sea_fin$lnonviolent - (mean(sea_fin$lnonviolent)/sd(sea_fin$lnonviolent))
sea_fin$lviolent_sd <- sea_fin$lviolent - (mean(sea_fin$lviolent)/sd(sea_fin$lviolent))
sea_fin$lsociety_sd <- sea_fin$lsociety - (mean(sea_fin$lsociety)/sd(sea_fin$lsociety))
sea_fin$lperson_sd <- sea_fin$lperson - (mean(sea_fin$lperson)/sd(sea_fin$lperson))
sea_fin$lproperty_sd <- sea_fin$lproperty - (mean(sea_fin$lproperty)/sd(sea_fin$lproperty))

# rename racial and economic boundary variable
sea_fin$sea_white_blv <- sea_fin$p_race_white_blv
sea_fin$sea_ses_blv <- sea_fin$ses_blv


#### reg tables associated with boundaryXp_white analysis for misdemeanor arrests ####

atl.race.int.misd = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                                age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                lviolentcrime + phys_bound + com_density, data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus.race.int.misd = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                                age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                lviolentcrime + phys_bound + com_density, data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos.race.int.misd = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                                age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, 
                              data = bos_fin, 
                              cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.race.int.misd = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                                age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                lviolentcrime + phys_bound + com_density, data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil.race.int.misd = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                                age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                lviolentcrime + phys_bound + com_density, data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea.race.int.misd = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                                age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                                lviolentcrime + phys_bound + com_density, data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# create latex table for boundaryXp_white for misdemeanors
texreg(l = list(atl.race.int.misd, aus.race.int.misd, bos.race.int.misd, chi.race.int.misd, mil.race.int.misd,
                sea.race.int.misd),
       file = "regtable.hetefx.misd.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "p_race_white" = "% White",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "lcrime" = "Log(Crime)",
                              "boundary_quart_dummy:p_race_white" = "Boundary X % White",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Percent White on Logged Misdemeanor Arrests (Standardized), Conditional on Racial Boundary Status",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("ATL","AUS","BOS","CHI","MIL","SEA"),
       custom.header = list("Dependent Variable: Misdemeanor Arrests" = 1:6),
       label = "table:misd_whiteXb",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))

#### reg tables associated with boundaryX\%white analysis for felony arrests ####

atl.race.int.fel = lm_robust(lfelonies_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                               age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                               lviolentcrime + phys_bound + com_density, data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus.race.int.fel = lm_robust(lfelonies_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                               age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                               lviolentcrime + phys_bound + com_density, data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos.race.int.fel = lm_robust(lfelonies_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                               age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lcrime + phys_bound + com_density, data = bos_fin, 
                             cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.race.int.fel = lm_robust(lfelonies_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                               age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                               lviolentcrime + phys_bound + com_density, data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil.race.int.fel = lm_robust(lfelonies_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                               age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                               lviolentcrime + phys_bound + com_density, data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea.race.int.fel = lm_robust(lfelonies_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + 
                               age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + lpropertycrime + 
                               lviolentcrime + phys_bound + com_density, data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# create latex table for boundaryXp_white for felonies
texreg(l = list(atl.race.int.fel, aus.race.int.fel, bos.race.int.fel, chi.race.int.fel, mil.race.int.fel,
                sea.race.int.fel),
       file = "regtable.hetefx.fel.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "p_race_white" = "% White",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "lcrime" = "Log(Crime)",
                              "boundary_quart_dummy:p_race_white" = "Boundary X % White",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Percent White on Logged Felony Arrests (Standardized), Conditional on Racial Boundary Status",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("ATL","AUS","BOS","CHI","MIL","SEA"),
       custom.header = list("Dependent Variable: Felony Arrests" = 1:6),
       label = "table:fel_whiteXb",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))


#### raceXclass plots ####
# now run interaction models with binned racial boundary value

### Atlanta ###
form_list = 
  list(
    'lmisdemeanors_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + age_15_35_male  + 
    pown  + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density',
    'lfelonies_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + age_15_35_male  + 
    pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density'
  )

hetdf_list = list(form_list)
hetdf_lab = c("Misdemeanors","Felonies")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), atl_fin)
  
  fakedata = 
    atl_fin[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Atlanta")

# create interaction plots
library(ggh4x)

atl_raceXclass = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_poverty, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% Poverty",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") + 
  theme_bw(base_size=12,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="none",
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.title.x=element_blank())



### Austin ### 

# now run interaction models with binned racial boundary value
form_list = 
  list(
    'lmisdemeanors_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + 
    age_15_35_male  + pown  + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density',
    'lfelonies_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + age_15_35_male  + 
    pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density'
  )

hetdf_list = list(form_list)
hetdf_lab = c( "Misdemeanors","Felonies")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), aus_fin)
  
  fakedata = 
    aus_fin[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Austin")

# create interaction plots

aus_raceXclass = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_poverty, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% Poverty",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") + 
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="none", # remove legend
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + # remove grey from background
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # remove gridlines from background
        axis.title.x=element_blank(), strip.text.x = element_blank()) # remove x axis title and x labels for plots


### Boston ###

# now create binned racial blv measures
bos_fin <- bos_fin %>% mutate(boundary_quart = quantile(bos_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))


# now run interaction models with binned racial boundary value
form_list = 
  list(
    'lmisdemeanors_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + age_15_35_male  + 
    pown  + p_emp_unemployed + pcol + lcrime + phys_bound + com_density',
    'lfelonies_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + age_15_35_male  + 
    pown + p_emp_unemployed + pcol + lcrime + phys_bound + com_density'
  )

hetdf_list = list(form_list)
hetdf_lab = c("Misdemeanors","Felonies")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), bos_fin)
  
  fakedata = 
    bos_fin[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Boston")

# create interaction plots

bos_raceXclass = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_poverty, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% Poverty",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x")  +  
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="none", # remove legend
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + # remove grey from background
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # remove gridlines from background
        axis.title.x=element_blank(), strip.text.x = element_blank()) # remove x axis title and x labels for plots




### Chicago ###

# now run interaction models with binned racial boundary value
form_list = list(
  'lmisdemeanors_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + 
  age_15_35_male + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density',
  'lfelonies_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + age_15_35_male + 
  pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density'
)

hetdf_list = list(form_list)
hetdf_lab = c("Misdemeanors","Felonies")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), chi_fin)
  
  fakedata = 
    chi_fin[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Chicago")

# create interaction plots

chi_raceXclass = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_poverty, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% Poverty",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") + 
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="none", # remove legend
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + # remove grey from background
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # remove gridlines from background
        axis.title.x=element_blank(), strip.text.x = element_blank()) # remove x axis title and x labels for plots


### Milwaukee ###

# now run interaction models with binned racial boundary value
form_list = 
  list(
    'lmisdemeanors_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + 
    age_15_35_male  + pown  + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density',
    'lfelonies_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + 
    age_15_35_male + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density'
  )

hetdf_list = list(form_list)
hetdf_lab = c("Misdemeanors","Felonies")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), mil_fin)
  
  fakedata = 
    mil_fin[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Milwaukee")

# create interaction plots

mil_raceXclass = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_poverty, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% Poverty",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x") + 
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="none", # remove legend
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + # remove grey from background
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # remove gridlines from background
        axis.title.x=element_blank(), strip.text.x = element_blank()) # remove x axis title and x labels for plots



### Seattle ###

# now run interaction models with binned racial boundary value
form_list = 
  list(
    'lmisdemeanors_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + 
    age_15_35_male  + pown  + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density',
    'lfelonies_sd ~ boundary_quart_dummy * p_poverty + p_race_white + ses_blv + lpop + hhi + lmhhi + 
    age_15_35_male  + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density'
  )

hetdf_list = list(form_list)
hetdf_lab = c("Misdemeanors","Felonies")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), sea_fin)
  
  fakedata = 
    sea_fin[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_poverty = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

df_to_plot = df_to_plot %>%
  mutate(city = "Seattle")

# create interaction plots

sea_raceXclass = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_poverty, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_poverty,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% Poverty",
       y = "Predicted Value",
       col = "Race Boundary") +
  facet_grid2(city ~ dataset, scales = "free_x", independent = "x")  + 
  theme_bw(base_size=12,base_family="Times") +
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom", # remove legend
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + # remove grey from background
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # remove gridlines from background
        strip.text.x = element_blank()) # remove x axis title and x labels for plots




# arrange all plots on one page (expect Louisville)

library(gridExtra)
library(cowplot)
all_cities_class_mod_plots <- plot_grid(atl_raceXclass, aus_raceXclass, bos_raceXclass,
                                        chi_raceXclass,mil_raceXclass,sea_raceXclass, 
                                        align = "v", nrow = 6, rel_heights = c(1/5, 1/6, 1/6, 1/6, 1/6, 1/5))

ggsave(plot = all_cities_class_mod_plots, filename = "all_cities_class_mod_plots.png", 
       width = 8, height = 14, bg = 'white')


## RaceXClass Analysis ##

#### reg tables associated with boundaryXp_poverty analysis for misdemeanor arrests ####

atl.misd.class.int = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
                                 hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                 phys_bound + com_density, 
                               data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)


aus.misd.class.int = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
                                 hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                 phys_bound + com_density,
                               data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos.misd.class.int = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
                                 hhi + lmhhi + pown + p_emp_unemployed + pcol + lcrime +
                                 phys_bound + com_density, 
                               data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.misd.class.int = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
                                 hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                 phys_bound + com_density, 
                               data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# lou.bound.class = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
#                               hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
#                             data = lou_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil.misd.class.int = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
                                 hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
                               data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea.misd.class.int = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
                                 hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
                               data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# create latex table for boundaryXp_poverty for misd arrests
texreg(l = list(atl.misd.class.int, aus.misd.class.int, bos.misd.class.int, chi.misd.class.int,
                mil.misd.class.int, sea.misd.class.int),
       file = "regtable.raceXclass.misd.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "p_race_white" = "% White",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "lcrime" = "Log(Crime)",
                              "boundary_quart_dummy:p_poverty" = "Boundary X % Poverty",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Percent Poverty on Logged Misdemeanor Arrests (Standardized), Conditional on Racial Boundary Status",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("ATL","AUS","BOS","CHI","MIL","SEA"),
       custom.header = list("Dependent Variable: Misdemeanor Arrests" = 1:6),
       label = "table:raceXclass_misd",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))


#### reg tables associated with boundaryXp_poverty analysis for felony arrests ####

atl.fel.class.int = lm_robust(lfelonies_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                phys_bound + com_density, 
                              data = atl_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

aus.fel.class.int = lm_robust(lfelonies_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                phys_bound + com_density,
                              data = aus_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

bos.fel.class.int = lm_robust(lfelonies_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + 
                                age_15_35_male + hhi + lmhhi + pown + p_emp_unemployed + pcol + lcrime +
                                phys_bound + com_density, 
                              data = bos_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

chi.fel.class.int = lm_robust(lfelonies_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                phys_bound + com_density, 
                              data = chi_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# lou.bound.class = lm_robust(lfelonies_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
#                               hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime, 
#                             data = lou_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

mil.fel.class.int = lm_robust(lfelonies_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                phys_bound + com_density, 
                              data = mil_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

sea.fel.class.int = lm_robust(lfelonies_sd ~ boundary_quart_dummy*p_poverty +ses_blv + lpop + p_race_white + age_15_35_male + 
                                hhi + lmhhi + pown + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime +
                                phys_bound + com_density, 
                              data = sea_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# create latex table for boundaryXp_poverty for misd arrests
texreg(l = list(atl.fel.class.int, aus.fel.class.int, bos.fel.class.int, chi.fel.class.int,
                mil.fel.class.int, sea.fel.class.int),
       file = "regtable.raceXclass.fel.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "p_race_white" = "% White",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "p_poverty" = "% Poverty",
                              "pown" = "% Homeowner",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)",
                              "lcrime" = "Log(Crime)",
                              "boundary_quart_dummy:p_poverty" = "Boundary X % Poverty",
                              "phys_bound" = "Physical Boundary",
                              "com_density" = "Commercial Density"),
       caption = "Influence of Percent Poverty on Logged Felony Arrests (Standardized), Conditional on Racial Boundary Status",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("ATL","AUS","BOS","CHI","MIL","SEA"),
       custom.header = list("Dependent Variable: Felony Arrests" = 1:6),
       label = "table:raceXclass_fel",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))

# For Louisville:
# Louisville: Influence of Percent White on Logged Arrests by Type (Standardized)
# now run interaction models with binned racial boundary value
form_list = 
  list(
    'lsociety_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + age_15_35_male  + pown + 
    p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density',
    'lperson_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + age_15_35_male  + pown + 
    p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density',
    'lproperty_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + lmhhi + age_15_35_male  + pown + 
    p_poverty + p_emp_unemployed + pcol + lpropertycrime + lviolentcrime + phys_bound + com_density'
  )

hetdf_list = list(form_list)
hetdf_lab = c("Against Society",
              "Against Persons","Against Property")

out_mod_list = as.list(rep(NA, length(form_list)))
out_df_list = as.list(rep(NA, length(form_list)))

for (i in 1:length(form_list)) {
  
  print(paste0("Iteration ", i))
  
  out_mod_list[[i]] = lm_robust(as.formula(form_list[[i]]), lou_fin)
  
  fakedata = 
    lou_fin[, names(out_mod_list[[i]]$coefficients)[2:(length(names(out_mod_list[[i]]$coefficients))-1)]] %>% 
    filter(lpop > 0) %>% 
    as.data.frame %>% 
    mutate(geometry = NULL) %>% 
    apply(X = ., MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE)) %>% 
    t %>% 
    as.data.frame %>% 
    slice(rep(1:n(), each = 22)) %>% 
    mutate(p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
           boundary_quart_dummy = c(rep(0, 11), rep(1, 11)))
  
  out_df_list[[i]] = data.frame(
    
    est = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$fit,
    se = predict(object = out_mod_list[[i]], newdata = fakedata, se.fit = TRUE)$se.fit,
    p_race_white = rep(seq(from = 0, to = 1, by = .1), 2),
    boundary_quart_dummy = c(rep(0, 11), rep(1, 11))
    
  ) %>% 
    mutate(dataset = hetdf_lab[i])
  
}

df_to_plot = out_df_list %>%
  do.call(rbind.data.frame, .) %>%
  mutate(boundary_quart_dummy = factor(boundary_quart_dummy, labels = c("No Boundary", "Boundary")))

# create interaction plots

lou_racial_blv_int = df_to_plot %>%
  ggplot() +
  geom_point(aes(x = p_race_white, y = est, col = boundary_quart_dummy),
             position = position_dodge(.075),
             size = 2) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .4,
                position = position_dodge(.075)) +
  geom_errorbar(aes(x = p_race_white,
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = boundary_quart_dummy),
                width = 0,
                size = .6,
                position = position_dodge(.075)) +
  scale_color_grey(start = .6, end = 0) +
  labs(x = "% White",
       y = "Predicted Value",
       col = "Race Boundary",
       title = "Louisville Moderation Analysis") +
  facet_wrap(~ dataset, scales = "free") + theme_bw(base_size=12,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA)) +
  theme(strip.background = element_rect(fill="white")) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())


ggsave(width = 8, height = 6, plot = lou_racial_blv_int, filename = "lou_racial_blv_int.png")


#### reg tables associated with boundary X %white --> arrests by UCR type Louisville ####

# Louisville
lou.person = lm_robust(lperson_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + hhi + 
                         lmhhi + age_15_35_male  + pown + p_poverty + p_emp_unemployed + pcol + 
                         lpropertycrime + lviolentcrime + phys_bound + com_density, 
                       data = lou_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

lou.property = lm_robust(lproperty_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + 
                           hhi + lmhhi + age_15_35_male  + pown + p_poverty + p_emp_unemployed + 
                           pcol + lpropertycrime + lviolentcrime + phys_bound + com_density, 
                         data = lou_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

lou.society = lm_robust(lsociety_sd ~ boundary_quart_dummy * p_race_white + ses_blv + lpop + 
                          hhi + lmhhi + age_15_35_male  + pown + p_poverty + p_emp_unemployed + 
                          pcol + lpropertycrime + lviolentcrime + phys_bound + com_density, 
                        data = lou_fin, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

# create latex table for boundary X %white --> ucr types of arrests for Louisville
texreg(l = list(lou.person, lou.property, lou.society),
       file = "regtable.hetefx.lou.tex",
       include.ci = FALSE,
       custom.coef.map = list("boundary_quart_dummy" = "Boundary (White)",
                              "p_race_white" = "% White",
                              "ses_blv" = "Boundary (SES)",
                              "lpop" = "Log(Population)",
                              "age_15_35_male" = "Age 15-35 Male",
                              "hhi" = "Diversity",
                              "lmhhi" = "Log(MHHI)",
                              "pown" = "% Homeowner",
                              "p_poverty" = "% Poverty",
                              "p_emp_unemployed" = "% Unemployed",
                              "pcol" = "% College",
                              "boundary_quart_dummy:p_race_white" = "Boundary * % White",
                              "lpropertycrime" = "Log(Property Crime)",
                              "lviolentcrime" = "Log(Violent Crime)"),
       caption = "Louisville: Influence of Percent White on Logged Arrests by Type (Standardized), Conditional on Racial Boundary Status.",
       include.rmse = FALSE,
       include.adjrs = FALSE,
       caption.above = TRUE,
       float.pos = "H",
       custom.model.names = c("Against Persons","Against Property","Against Society"),
       label = "table:lou_whiteXb",
       symbol = "\\dagger",
       stars = c(0.001, 0.01, 0.05, 0.1))



